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Abstract 

We obtain systematic approximations for the modes of vibration of a string of 
variable density, which is held fixed at its ends. These approximations are obtained 
iteratively applying three theorems which are proved in the paper and which hold 
regardless of the inhomogeneity of the string. Working on specific examples we 
obtain very accurate approximations which are compared both with the results of 
WKB method and with the numerical results obtained with a collocation approach. 
Finally, we show that the asymptotic behaviour of the energies of the string obtained 
with perturbation theory, worked to second order in the inhomogeinities, agrees with 
that obtained with the WKB method and implies a different functional dependence 
on the density that in two and higher dimensions. 

Key words: Helmholtz equation; inhomogeneous string; perturbation theory; 
collocation method 



1 Introduction 



We consider the problem of describing the vibrations of a string of variable 
density, which is held fixed at its ends (Dirichlet boundary conditions). This 
problem has been investigated in depth in a series of papers by Krein, who 
has considered both the direct and inverse problem [T|2|3]IHl5] . In recent years, 
Beals, Sattinger and Szmigielsky [6], have studied the string density problem 
in connection with the Camassa-Holmes equation for shallow water waves. 

Although the problem for a string is considerably simpler than the correspond- 
ing two dimensional version, i.e. an inhomogeneous membrane of arbitrary 
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density, exact solutions are known only in few cases: two examples of exactly 
solvable problems are the strings studied by Lord Rayleigh [7] and in the late 
40s by Borg [8J, which correspond to case where the reciprocal of the densi- 
ties vary as the fourth or the second power of the position respectively. Other 
examples of exactly solvable inhomogeneous strings are given by Horgan and 
Chang ig. 

Recently Gottlieb [10] has studied these examples and he has obtained a trans- 
formation which maps the original variable density string into a family of 
strings which are isospectral to it: in this way Gottlieb shows that the case 
discussed by Borg can be obtained from a homogeneous string via his trans- 
formation. Another case of problem which can be formally solved exactly is a 
string of density which depends linearly on the position: Fulcher pJJ has ob- 
tained the exact solutions for this case, in terms of a complicated trascendental 
equation. 

In the cases which cannot be solved exactly or where the exact solution is of 
no practical use, as for the example discussed by Fulcher, one needs to resort 
to approximations: a typical approach is the WKB method, which correctly 
describes the leading behaviour of the highest part of the spectrum, although 
it is less precise for the lowest excited modes 0- Using this approach, Crawford 
has obtained a simple analytical formula for the problem solved by Fulcher, 
which describes all the spectrum with an accuracy of few percent [13]. An 
example of application of the WKB method beyond the leading order to a 
nontrivial problem of a density with a rapidly oscillatory behaviour is dis- 
cussed by Castro and Zuazua [TifTS] . An interesting computational scheme 
for the solution of the inhomgeneous one dimensional Helmholtz equation has 
been recently discussed by Rawitscher and Liss in ref. [15], using a spectral 
expansion in term of Chebyshev polynomials. 

An alternative to the WKB method is provided by perturbation theory (PT), 
which however is bounded to cases where the density is only slightly perturbed 
with respect to the density of an exactly solvable problem. Using this approach 
in the case of two dimensional membranes (or billiards), and performing a 
suitable resummation of the perturbative terms, it has been possible to derive 
Weyl's law, which relates the spectrum of the billiard to its area (see [IZ]). In 
the present paper we show an interesting (and unexpected) result: the Weyl's 
law for a one dimensional string corresponds to a different functional of the 
density than in higher dimensions (the problem of a homogeneous membrane in 
two dimensions is isospectral to the problem of an inhomogeneous membrane 
with density obtained from the conformal map which send the border of the 
membrane to a reference border). The results that we obtain here, which are 



The application of the WKB method to the variable density string is discussed 
in detail by Bender and Orszag in their book [12j . 
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exact to second order, agree with the general WKB formula and provide a link 
between these two different approaches. 

The central results of the present paper are contained in three theorems, which 
provide an explicit iterative tool to build increasingly accurate approximations 
to the modes of the strings, without being restricted to slightly inhomogeneous 
systems or highly excited states. We discuss the application of these theorems 
to specific examples. The comparison of these results with the precise numeri- 
cal results obtained with a collocation method allows to verify explicitly their 
convergence. 

The paper is organized as follows: in Section |2] we describe the collocation 
approach; in Section [3] we enunciate and prove the theorems; in Section H] 
we discuss the application of Perturbation Theory and prove that, to second 
order, one recovers the results of WKB method; in Section [5] we discuss two 
applications of the general results obtained in the previous sections; finally, in 
Section [H] we draw our conclusions. 



2 Collocation approach 

In this section we briefly describe a collocation approach which can be used 
to solve the Helmholtz equation for a one dimensional string with variable 
density. The method that we are using is the Conformal collocation method 
(CCM) that we have devised in Ref. [IB], where we have shown that it can 
be used to obtain precise numerical approximations to the energies and wave 
functions of the Helmholtz equation on an arbitrary two dimensional region 
with Dirichlet boundary conditions. This method has also been used recently 
to obtain very precise numerical estimates for arbitrary domains in the plane, 
either simply connected, P^|IT7] . or with a hole, [2Dj, these last domains being 
known as " quantum rings" . 

In order to make our discussion self contained we briefly review the main 
features of the method. Our starting point is the inhomogeneous Helmholtz 
equation on the interval [—L,L]: 



where p{x) > is the density of the string and ■?/'„(±L) = 0. The eigenfunctions 
ipnix) are orthogonal with respect to the weight function 




(1) 
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L 

J llJn{x)lpra{,x)p{x)dx = 5nm ■ (2) 
-L 



In order to discretize this eigenvalue equation we introduce a set of functions, 
which we call Little Sine Functions (LSF), discussed in [21] and fuUfilling 
Dirichlet boundary conditions 0: 



N sin (if ) - sin (f ) • 



with k = —N/2 + 1, . . . , N/2 — 1. These functions define an homogeneous grid 
2Lk 

Xk = — ^ and obey the orthogonality relation 



Sk{N, L, x)sj{N, L, x)dx = h6kj, (4) 



-L 



where h = ^ is the spacing of this grid. 



N 



A function f{x) obeying Dirichlet boundary conditions may be interpolated 
using the Sk{h, N, x) as 



N/2-1 

fix)^ E f{xk)skih,N,x) . (5) 

k=-N/2+l 



To better understand the validity of this expression we recall the definition of 
LSF 



2L ^ 

Sk{N, L,X) = 1pn{x)iJniXk) , (6) 

n=l 

from which eq. was obtained in ref.[2I]. Here ipn{x) = ^ sin "^(^+^) are 
the solutions for a homogeneous string of unit density, which form a basis. 

Substituting this expression inside the rhs of eq. ([5]), we obtain 



^ In [_22j we have also discussed three other sets of LSF obeying different boundary 
conditions: choosing one of these alternative sets would then correspond to solve 
Helmholtz equation ^ with the corresponding boundary condition. 
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N 

E 

n=l 



2L ^''^"^ 
k-N/2+1 



(7) 



One may easily recognize that the term in parenthesis is a Riemann sum and 
that, for — )■ oo, it converts to the integral J^l f{y)ipn{y)dy, thus leading 
to the standard decomposition of a function f{x) in the basis (x)}. This 
discussion can also be found in Ref. [2T]. 



In the same way one may obtain an interpolation formula for the second 
derivative of the function, simply deriving twice the expression above and 
thus obtain 



X) 



dx'^ 



N/2-l 



<Psk{x) 
dx"^ 



d^Sk{x) 



dx"^ 



k=-N/2+l 

N/2-1 N/2-1 

' E E 

fc=-Ar/2+l j=-Af/2+l 
N/2-\ N/2-1 

--- E E /(x.)cg).,(/i,iv,x), 

fc=-Ar/2+l j=-7V/2+l 



Sjih, N, x] 



where the matrix of elements c^.^"* = ^^i^"^ 
second derivative operator on the grid. 

After writing eq. (JT]) into the equivalent form 



provides a representation for the 



1 d^ 

p{x) dx' 



:ip{x) = Ei/j{x) 



(9) 



one may obtain a representation of the differential operator on the left, O = 
— 7-T-7-2, in terms of the LSF introduced above as: for a set of LSF with a 

given N, the matrix elements of O on the grid are given by 



Okk' 



.(2) 



Pixk^'"' 

where k, k' = -N/2 + 1, . . . , N/2 - 1. 



(10) 



As pointed out in [19] this form of O is not manifestly hermitean, so that 
it is more convenient to work with the symmetrized form of the operator 



O 



sym 



^ p{x) "'^^ a/ 



whose matrix elements read 
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n'^y"' - p(^) nil 



The only differences in this discussion from the cases discussed in 
are that we are hmiting ourselves to one dimension and that p{x) is a physical 
density, i.e. cannot be considered as obtained from a conformal mapping. 

The reader may be worried that dealing with the symmetrized operator may 
not be equivalent to dealing with the original operator: we may easily convince 
ourselves that this is not so. Let us consider the Helmholtz equation iQ and 

introduce the function (f){x) = ^ p{x)ilj{x). If we substitute this function in 
the equation we obtain 



[x) = E<P{x) , (12) 



which involves the symmetrized operator introduced earlier. In essence the 
eigenvalues of O and Ogym are the same, while the eigenfunctions are related 



by a factor p{x): this also implies that the eigenfunctions of eq.(|9]) are or- 
thogonal with respect to the weight p{x), as stated in eq. (I2l). 

The advantages of the present collocation approach are clear: 

• the representation on the grid of the Helmholtz equation with variable den- 
sity is obtained straightforwardly; 

• the matrix representing the differential operator is the product of a universal 
matrix (the matrix for the laplacian on a finite interval) with the matrix (or 
matrices depending if one is using the symmetrized form of the operator or 
not) which depends on the density but is diagonal; 

• the calculation of the universal matrix corresponding to a given grid may 
be done once and for all, whereas the diagonal density matrices need to 
be calculated for each specific density; clearly the computational price of 
this second operation is negligible with respect to the first task. The correct 
computational strategy can therefore be calculating the matrix of the Id 
laplacian for different grids and then store them for later use; 

• at no time has one to perform numerical integrations; 

• the collocation method provides exact results for constant densities: for this 
reason its application to problems with slightly inhomogeneous densities 
provides very precise results; 

After this process has been carried out one obtains a {N — 1) x {N — 1) her- 
mitean matrix whose eigenvalues and eigenvectors will provide approximations 
to the lowest — 1 energies and wave functions of eq. (P). 
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3 Non-perturbative approach 



In this section we wish to obtain explicit formulas for the states of a string of 
variable density which obey eq. ([1]). These formulas are non-perturbative, i.e. 
they do not depend on the inhomogeneities being small. 

We state the following theorem: 

Theorem 1 Let S,q{x) he an arbitrary function defined on the interval {—L, L) 
(being C,o{x) arbitrary we assume that it has nonzero overlap with the true 
fundamental solution). For n — )■ oo the sequence of functions 



p{x) J dy 

-L 



dz\[p{z)in-l{z) 



-L 



(13) 



with K„ = 2X S-L^y I-L dz \J p{z)^n-i{z) , convcrgcs to the fundamental mode 
of the Helmholtz equation (Qp with Dirichlet boundary conditions, ^^(ibL) = 0.' 



*o(a;) = lim in{x 

n^oa 



(14) 



where Cn{x) = ■^M=, 



Proof Since ^0(2;) is an arbitrary function of x, it will have a nonzero overlap 
with the exact wave function of the ground state, '^q{x). The operator O = 
— -^-^-^ associated to the eq. has a spectrum which is bounded from 

below, positive and simple (Theorem 2.1 of [6]). The spectrum is not bounded 
from above. The inverse operator therefore has a spectrum which is also 
positive and simple, but it is now bounded both from below and from above. 
The function obtained by applying to ^o{x), corresponding to the case 
n = 1 in eq. (fT3|) obeys Dirichlet boundary conditions and it has a larger 
overlap with the exact wave function ^^0(2;) O 



J^l: ./^^oix)Ux)dx\ ^ ,/^^o{x)Ux)dx 



(15) 



As a matter of fact the operator O ^ acting on ^^{x) "inflates" the compo- 
nents of ^q[x) corresponding to the lowest energy states of O. The iteration 



Recall the orthogonality relation ([2]). 
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of this precedure will then provide a sequence of functions which converge to 



'^q{x)\ p{x) as n — )■ oo. In the case that io{x) is orthogonal to the first k 



states of eq. ([I]), ^q{x)^J p{x), . . . , \E'fc_i(a;)yp(a;), then the sequence of func- 
tions ^„(x) will converge to ^ k{.x) \J p{x) . 



Incidentally, one may consider this theorem (and the more general Theorem 1 
of [19]) as the generalization of the celebrated power method of linear algebra 
to operators in a Hilbert space. Notice also that eq. (fT3|) . where ^n{x) = 
= ^(a;), corresponds to Helmholtz equation in its integral form. 

We may easily extend this theorem to calculate the first states of the inho- 
mogeneous Helmholtz equation: 



Theorem 2 Let 2^°) = [s}^ 



X) 



,S,Q^\x)^ he a set of arbitrary functions 



defined on the interval {—L, L). We assume that these functions have a nonzero 
overlap with the first N solutions of the inhomogeneous Helmholtz equation. 
For convenience we assume that the functions are ordered as 



l!:Ld'\x)'dx 



I-L6''\^)0eo"'ix)dx 



< 



xYdx 



(16) 



Consider the sequence of functions 



^^\x) ^ VP(^ j dy 

-L 



4^'^ - j dz^p{z)e^\z) 

-L 



With kY' = 2L I-idy I-Ldz^p{z)Co {z) and] = 1, 
Let 



N. 



(17) 



j-l fL Ak)(\A3) 



i-Lcrim"{x)dx- 



k=l 



f':^if\xydx 



and S*^^) = ^^i'\x) , . . . , ^i^\x)^ . We call 2*^"^ the set of functions obtained 

after n iterations. Then, for n — )■ oo, S*^") converges to the N lowest eigen- 
f unctions of eq. (QP- 



Proof The proof is straighforward: at each iteration, the functions generated 
with (ITTl) have their lowest energy components inflated. The orthogonalization 
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in eq. ( !T8|) eliminates from the j*'* functions the components corresponding to 
the previous j — 1 functions, which have lower expectation values. As this 
procedure is iterated only the j*'^ exact mode will survive in the j*'* function. 
This completes the proof. 



The two theorems that we have just proved allow to build iteratively the 
solutions corresponding to the lowest modes of the inhomogeneous string. 
A different strategy to calculate the ground state of the string consists of 
using eq. f[T^ in a non iterative way applying the Variational Theorem: if we 
let ^o{x) be an arbitrary function, depending on one or more parameters Ui, 
then one may determine these parameters minimizing the Rayleigh quotient 
corresponding to C,i{x), i.e. the function obtained with eq. ( !T3|) applied to 
^o{x)- Alternatively one could minimize the Rayleigh quotient corresponding 
to C,o{x), but this is clearly less efficient, since Ci(^) is always closer to the 
fundamental mode because of Theorem [TJ 

We may also formulate a non-iterative procedure which allows to obtain ap- 
proximations for the lowest part of the spectrum of the string. We now briefly 
describe this approach. 

Let H'^^) = ^^o^\x), . . . ,^o^\x)^ be a set of functions defined on the interval 
(— L, L), and let S*^^^ = \^^i'\x), . . . , ^i^\x)^ be the set of functions obtained 

after applying eq. (fT3|l to each of the functions in and normalizing each 
of them: 



Ci\xfdx = l. (19) 



Let A and B be the matrices whose elements are 

^ ( 1 ^ eP^(x) (20) 

B.,= / ■ (21) 




The eigenvalues A and eigenvectors v of the generalized eigenvalue problem 
Av = ABv (22) 
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are approximations to the lowest modes of the string. The advantage of this 
non-iterative procedure is that the set Sq may be chosen arbitrarily (for ex- 
ample, selecting functions for which the integrals can be done analytically) 
and that there is no need of orthogonalization, as in Theorem |2J 

A yet different approach would consist of applying this procedure to a set of 
functions Sq = {1^0(2^)5 • • • j^n{x)}, obtained after N iterations of eq. flT^ on 
an arbitrary function ^0(2;) (we also assume that each function is normalized). 
In this case So(x) defines a Krylov subspace, and it can be used to set up a 
generalized eigenvalue problem as in eq. ( 122|) . One expects this procedure to 
be more efficient that the previous, given that the iteration of eq. f ll3p pro- 
vides functions whose lower energy components are enhanced: however, from a 
practical point of view, the iteration of eq. f lT3|) in general generates functions 
which are more and more complicated (this problem could be obviated, at 
least in part, by finding suitable analytical approximations for each function). 

The results that we have described so far allow to obtain arbitrarily precise 
approximations for the lowest modes of an inhomogeneous string. 

We now enunciate a third theorem, which allows one to calculate excited states 
of the string without the need of orthogonalization with respect to the lower 
lying modes. 

Theorem 3 Let riQ{x) be an arbitrary function fulfilling Dirichlet boundary 
conditions on [—L,L'], with — L < L' < L. The sequence of functions 



Vnix) = yp(^ / dy 



dz^J p{z)rin-i{z) 



(23) 



with Kn = x+l^ I-L I-L ^'^\J p{z)Vn-i {z) , convcrgcs to a solution of equation 
(Cp for n ^ 00 if L' is chosen so that rjn{L) = 0. 



Proof The proof goes as follows: for a fixed L' the sequence of functions 
r]n{x) would converge to the ground state of the Helmholtz equation ([1]) with 
Dirichlet boundary conditions in x = —L and x = V because of Theorem [H 
By fixing V at each iteration so that rin{L) = 0, one obtains a function which 
converges to a solution of eq. ([T]) on the whole string. 



Notice that as the number n of iterations is increased, the equation rin{L) = 
acquires multiple solutions, each one approximating a different solution of 
eq. (P). 
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We will apply these theorems to specific examples in Section |5l 



4 Perturbation theory 



While in the previous sections we made no assumption regarding the behaviour 
of p{x), apart from requiring it to be a positive function, now we also assume 
it to represent a slightly inhomogeneous string . In this case the problem 
can be treated with the perturbative approach of ref. [19]. Notice that the 
perturbation method developed in that paper may be regarded as a shape 
perturbation method, since the " density" was obtained from a conformal map 
and therefore changes in the density corresponded to changes of the shape 
of the membrane. In this case however p{x) is really a mass density, so that 
we may regard this specific application of perturbation theory as a genuine 
example of density perturbation theory. 

We briefly review here the perturbative approach described in ref. [TSl] and 
adapt it to the specific problem that we are considering: the starting point is 
the symmetrized operator 



= ^(-A)^, (24) 

which reduces to the negative laplacian for a constant p. In the present case 
it is understood that A is the second derivative (P/dx^. 

We now express p(x) as 



p{x)= pQ{l + 7]a{x)) (25) 

where \o{x)\ ^ 1 and po = ^I^LPi^)^^ is the average mass density in 
{—L,L). The parameter t] has been introduced for power counting; we will 
use T] to keep track of the different orders in powers of a{x) and set it to 1 at 
the end. 

Our operator may now be expanded in powers of rj as 
1 



Po 



Oo + vOi + v'02 + v'03 + ...\ , (26) 
where the explicit form of the may be worked out rather easily: 

Oo = -A (27) 
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6i = --[a(-A) + (-A)a] 



O2 
O3 



1 r 



2(t(-A)ct + 3(t\-A) + 3(-A)(T^ 



3 r 



16 



a\-A)a + a{-A)a' 



16 



a%-A) + {-A)a' 



(28) 
(29) 

(30) 



This problem may be treated within the standard Rayleigh-Schrodinger per- 
turbation theory (RSPT), keeping in mind that the operator O contains the 
perturbation to all orders in 77 (typical applications of RSPT see the pertur- 
bation term in the Hamiltonian to be of order one in the coupling constant). 

Working up to third order one finds 



Po Ej^^ = en 
po^i'^ = (ri|0i|r2) 

poi^f = (n|02|n) + E 



MO, 



p, E^~^ = {n\6,\n) + 



{n\d2\k){k\di\n) 

^k 



_l_ ^ ^ {n\0i\m){m\0i\k){k\0i\n) 



in 



Oi\n)Y: 



kf^n 



(n|Oi|A;)- 



(31) 
(32) 

(33) 



(34) 



where 



4L2 



TT and \n) are the eigenvalues and eigenstates of —cP/dx^ 



As shown in [19] one may work out the explicit form of these equations and 
obtain the formulas: 



Po 
Po 



= —en{n\a\n) 



Po Ei'^ = en{n\a\ny + elYl 



kf^n 



{n\a\k)^ 



(35) 
(36) 

(37) 



Po 



E^^^ = -en{n\a\nf + ei{n\a\n) ^ 



{n\a\ky 



kf^n 



nk 
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Q 2/ I I \ {n\(T\kf 



s (nlalk) (k\a\m) (m\a\n) .^^ 



Unlike in the cases studied in ref. [I9] these formulas apply to all the spectrum, 
which in one dimension is nondegenerate. 

If one neglects in these formulas the terms mixing different states, one may 
see that the remaining terms correspond to the terms of a geometric series 
and obtain the approximation (see |19j ) 



^ = ^" . (39) 

{n\po{l + (T{x))\n) {n\p{x)\n) ' 



This resummation was used in [T7] to rederive Weyl's law for two dimensional 
membranes and to generalize the law to two dimensional membranes of ar- 
bitrary density and shape (and in higher dimensions to cZ-cubes of arbirtrary 
density). We do not have a general proof that the terms mixing different states 
can always be neglected, although the fact that Weyl's law is obtained in this 
limit may support this claim. In the case of the small deformations of a square 
drum, we have proved in [19] that the contributions which mix different states 
at second order in perturbation theory are subleading with respect to the 
" diagonal" contributions. 

We will now show that in one dimension the Weyl's density law is modified. 
To start with we notice that for n — 00 



1 ' 



{n\a{x)\n) ^ — j a{x) dx . (40) 



-L 



Notice that in two dimensions, in the cases in which a{x, y) is related to a 
conformal density, one recovers in this limit the area of the membrane (see 
from which Weyl's law follows. 



Let us now focus of the second contribution in eq. (!37|) . which contains matrix 
elements of a between different states: 



E ■ (41) 
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We are interested in the behaviour of this term for n — i- oo. 
Let us consider an arbitrary cr(x) and write 



1 

{n\a\k) ~ ^ / cos 



[n — k)7r 
2L ' 



x + L) 



a{x)dx 



(42) 



which holds for n — > cxd. 
Therefore we have 



{n\a\k)'^ 
e„2^ ~ 



+L +L 

dx J dy a{x)A{x,y)a{y) 

L -L 



(43) 



where 



A{x,y) 



E 



n 



cos 



[n — k)'7i 
2L ' 



x + L) 



cos 



(n — k)7T 
2L ' 



X + L)(44) 



It is convenient to decompose cr(x) in a even and odd functions and therefore 
write 



{n\a\ky 



L +L 

dx / dy ae{x)A'^''\x,y)ae{y) + (Joix)A^°\x,y)ao\ 



where the sum in A'^'^\x,y) {A'^°\x,y) ) contains the values oi k ^ n, with 
\k — n\ even (odd). 

Let us first consider the function corresponding to odd indices: 



A^°\x,y)^Y. 

j=0 



cos cos 



(2j+l){L+y) 
2L 



8L2 



oo 

— 5]0j(x)(/)j(?/) 



(46) 



where (f)j{x) are the odd functions of the orthonormal basis fulfilling Neumann 
boundary conditions on the interval [—L,L] (see eq. (20) of [22]). Keeping in 
mind that A^°\x, y) only acts on odd functions, we may use the completeness 
of this set to make the following identification: 
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(47) 



We may now come to the function corresponding to even indices: 



A(^)(x,y)^E 
1 



oo COS (l^) COS 



8L2 



(48) 



where ipj{x) are the even functions of the orthonormal basis fulfilhng Neumann 
boundary conditions on the interval [— L, L] (see eq. (19) of [22] )• To exploit the 
completeness relation of this basis we must rewrite this expression including 
the zero mode as 



2L 



+ 



j=0 

6{x - y) 



16L2 ■ 8L 

keeping in mind that this function operates only on even functions. 



(49) 



We may use the results to write the asymptotic behaviour of the second order 
perturbative term: 



^(2) 



I \2 (n\a\k)'^ 
a\ny + enJ2 _ 

2 



— f (j{x)dx + / ci'^{x)dx— — f a{x)dx 
zL/ J 8-L J ALj J 



(50) 



Notice that this behaviour is compatible with the form of Weyl's law for a 
string of variable density p(x) = po(l + <^{x)) of the form 



E ^ 



p{x)y 



(51) 



where {y/p) = j^i \J p{x)dx. As a matter of fact we may expand eq. (!5T]) in 
powers of a and obtain: 
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7-1 II' 

En^ — 

Po 



l-{a) 



(52) 



which to second order reproduces the perturbative result that we have just 
discussed. Actually, one may derive this result using the WKB method, as done 
by Bender and Orszag in their classical book, [12], eq.(10.1.31) pag.490. For 
the case of a string of density which varies linearly with the position, Crawford 
[13] has compared the WKB formula with the exact results obtained by Fulcher 
for this particular problem, see [H], showing that it is very accurate. 



In a similar way one may also obtain the perturbative corrections to the states 
of the string: in this case, working to first order one obtains 



(0) (0) 



\k) 



_.(0) 



{k\(7\n) 
JO) JO) 



E 



(53) 



(54) 



where \n) are the unperturbed states. 



A general expression for the solutions of a string of variable density obtained 
with the WKB method is also given in eq. (10.1.33) of ref. [12] : 




-1/2 



-1/4 



sm 



niT- 



i-L \^pm 



(55) 



which is valid for n — )■ cxd. 



5 Applications 



In this Section we consider two examples of strings of variable density and 
apply to these problems the numerical and analytical results obtained in the 
previous sections. We adopt the convention of working with strings of unit 
length, centered in the origin. 
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5.1 "Borg string 



The first example that we consider has been originally discussed by Borg in 
[H] and corresponds to a density 



where a > —1 is a free parameter and 2;g( — 1/2,1/2). This string is isospec- 
tral to a string of uniform density p{x) = 1, and its solutions are known 
explicitly. 

The conditions for the isospectrality of two strings of different density have 
been recently discussed by Gottlieb in ref. [10], where the "Borg string" is one 
the examples considered in that paper. 

Gottlieb shows that the nonuniform Helmholtz's equations 



—i;(^) = Ep,iOm (57) 
Pix) = Ep2ix)(j){x) , (58) 



dx^ 

are isospectral if ^(x) = ^l^ai^l\~^2) ~^ \ ("^i^h a > — 1 and x,^ E (—1/2, 1/2)) 
and 



P2(x) = (f \ pi(e(x)). (59) 



The reader should notice that this relation is fully consistent with the asymp- 
totic law in eq. (ISTl) . since it implies 



1/2 1/2 



Ip2{x)dx= J ^pi(Oc?e (60) 

-1/2 -1/2 



In the case of a uniform string, the rational transformation considered by 
Gottlieb, provides precisely the density of the "Borg string". 

For this problem the asymptotic WKB method provides the exact energies 
and solutions; as a matter of fact, in the case of the energy, eq. (ISTll . reduces 
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to En = n^TT^, since (a/p) = 1; in the case of the wave functions, the WKB 
formula of eq. (!55|) reduces to 



/2{a + l) \ 2ax + a + 2 J 



which agrees with eq. (4.5) of Gottheb 2002 apart from a normahzation factor. 

We will now use the Borg string to test the numerical and analytical methods 
of our paper. 

We first focus on the application of the collocation method of Section [21 in the 
left plot of Figure [1] we display the error S„ = E^'^'^'^ /E'ff°''^^ — 1, for a = 1 (solid 
line) and for a = 10 (dashed line), where E^^'^^ are the energies calculated 
with the collocation method of section [2] using a grid corresponding to = 
2000. Notice that for a = 10 we are considering a strongly inhomogeneous 
string and the precision of the numerical results is clearly smaller, compared 
to the precision of the results corresponding to a = 1. An intuitive justification 
of this phenomenon is the following: as a grows, the density becomes strongly 
peaked around on of the ends of the string; in this case the excited modes 
of the string tend to oscillate more rapidly in this region and more slowly in 
the lower density region, compared to the corresponding modes of a uniform 
string. In the collocation approach, the number of grid points determines the 
maximal resolution which can be achieved in the calculation: for this reason 
the numerical precision obtained for the inhomogeneous modes, which need 
better resolution, is necessarily poorer. 

In the right plot of Figure [1] we compare the exact and numerical solu- 
tion for the fundamental mode of the Borg string for a = 10. The two 
curves are not distinguishable: the integrated error for the solutions is H = 
J-L 'ipi'"^\x) — ijjf^"''^^{x) ~ 3.36 X 10""^°, which is consistent with the error 
over the energy for the fundamental mode in the left plot of the figure. 



We would like to draw the attention of the reader on the flexibility of the 
collocation method: the calculation for strings of different density can be made 
quite efficiently specifying the density profile and evaluating it on an uniform 
grid. The construction of the diagonal density matrix is extremely fast and 
it takes typically a fraction of a second on an average desktop computer; the 
construction of the non-diagonal matrix for the laplacian is more involved, 
although it is general and it can be calculated once and stored. 

We now come to illustrate the application of Theorem [1] We pick the function 
^o{x) = 1 and apply to it eq. (|T3l) finding 
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Fig. 1. Left plot: Error over the energies of the Borg string, 



jpuuM I Tpexact 



1, 



^CCM 

for Q = 1 (sohd hne) and for a = 10 (dashed Une), calculated with a grid correspond- 
ing to = 2000. Right plot:Solution for the fundamental mode of the Borg string 

for a = 10. The dashed line is the exact result '4>i{x) = 2y^(5x+3) sin ^ ^4(^5^+3!'' ) • 
The solid line is the numerical result obtained with the collocation method using a 
grid with N = 2000. 



2{a + lf{{2x + 1) log(a + 1) - 2 log(2ax + a + 2) + log(4),) 



a^{2ax + a + 2)2 



which now vanishes at the ends of the string, ^i(±L) = (notice that this 
function is not normilized). We may easily iterate eq.f llSp and calculate explic- 
itly the next functions ^2(2;), ^3(2;), .... No numerical approximation is made 
in this process so that one does not need to worry about round-off errors. In 
this way we have obtained the expressions up to ^5(0;), although we do not 
report them here given their lengthy expressions. 

One may now consider the normalized function corresponding to the j^^ it- 
eration and decompose it in terms of the exact solutions, which form an or- 
thonormal and complete basis: 



k=l 



(63) 



where ^j{x) = ^j{x)/y J_]^^j{x)dx. The coefficients of this expansion are 



p{x) ^j{x)ipk{x)dx 



(64) 



with Ylh=i = 1 (Parseval relation). 
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3 


4 


5 


j 


0.59542214 


0.98958778 


0.99947247 


0.99996837 


0.99999804 


0.99999988 



Table 1 



First coefficient of the expansion of the functions ^j{x) corresponding to different 
iterations of eq. (|13p for the Borg string. 
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Fig. 2. Coefficients n'p of the expansion of the j*^ iteration $,j{x) in terms of the 
exact solutions of the Borg string for a = 10. 

Figure [2] and Table [1] illustrate the convergence of the {x) towards the fun- 
damental mode of the string: at each iteration the component of the function 
obtained from eq. (fT3l) corresponding to the fundamental mode is inflated with 
respect to remaining components and the approximate solution gets closer to 
the exact solution. 



We may calculate the expectation value of O = 
states corresponding to different iterations. For instance: 



in each of the 



di[aj + d2{a) log(l + a) + d^laj log (1 + a) 

— 1j 



where 



ni(a) = 108a^ , n2{a) = 108a\l + a) 

di{a) = (^a^ + 3a + 3^ , d2{a) = —30a (^a'^ + 3q; + 2 

4(0) = 36(1 + . 



(66) 
(67) 
(68) 
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2 

36.7836 


75117 
7550 

0.807197 


9800153 
992524 

0.0442409 


88566604437 
8973430934 

0.00269973 


1113581185772498961 
112829174318220286 

0.000167958 



Table 2 



Expectation values of the operator O in the states corresponding to different itera- 
tions and corresponding error. 

It is interesting to notice that this expression has a finite hmit both for a — )■ 
— 1 and for a — )■ 00, which correspond to the worse case scenarios for this 
problem. Notice that the two limit are clearly the same since the density is 
invariant under the transformations x — )■ —x and a — )■ —a/{l + a): therefore 
the Helmholtz equation ([1]) is also invariant under these transformations. The 
reader may explicitly check that the exact solutions of eq. (!6T]l are unchanged 
under the transformations above. 

As a matter of fact we have 



lim Ek'>(a) = lim Ek'>(a) = — , (69) 

which should be compared with the exact result E^^"-^^ = vr^. The first iteration 
therefore provides an estimate for the energy of the fundamental mode with 
an error which is at most of about 37% (for the case a = 10 studied before 
the error is of about 10%). 

The expressions corresponding to higher iterations can also be calculated ex- 
plicitly, although we do not report them here; in Tabled we display the ener- 
gies corresponding to a = — 1 and a = 00 for the different iterations and the 
corresponding error (in %) calculated with respect to the exact result. 

A different strategy for calculating the fundamental mode of a nonuniform 
string is to use the variational theorem. We pick an arbitrary function C,o{x) = 
1 + vx, where t; is a variational parameter which can be fixed by minimizing 
the expectation value of O in the state ^i(a;) obtained applying eq. (IT^ to 
^o{x)- In Fig. [3] we compare this variational energy with the energy obtained 
setting the parameter v = 0, which corresponds to the case studied earlier. The 
variational energy leads to a consistent gain in precision for moderate values 
of a, although for a — — 1 and a — 00 the two approximations provide the 
same limit, meaning that the term vx in ^o{x) is inessential in this regime. 
Clearly, this problem may be easily obviated by using different variational 
ansatz with the correct asymptotic behaviour and/or by using ansatzes with 
multiple variational parameters (possibly, simple enough so that the integrals 
in eq. (IT^ can be done explicitly). 

In Fig. m we illustrate Theorem [2] on the first 5 states of the Borg string, for 
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Fig. 3. Variational energy for the ground state obtained using £,o{x) = 1 + vx (solid 
curve) compared with the exact result (horizontal solid line) and with the result 
obtained setting to zero the variational parameter; the upper horizontal line is the 
upper bound 27/2. 



a = 10. We have used as starting functions the functions 



'5^ 



^Oy^)- I r(n+6) ' ^'^^ 

t+|)(n+l)! 



where C^^{{2x) are the Gegenbauer polynomials (we have chosen these func- 
tions because they are simple enough to perform explicit calculations). These 
functions are not a particularly good ansatz, as we can appreciate by looking 
at the open circles in the Figure, which show that each of these functions has 
an average energy which is far higher than the corresponding exact values. At 
each iteration however the ratio decreases and tends monotonically to 1 for all 
the states. From a practical point of view this Theorem can be used to obtain 
approximations for the few lowest energy modes, given that each new state 
needs to be orthogonalized with respect to the previous one, thus leading to 
a rapid increase in the number of operations needed. 



5.2 Parabolic density 



We consider now the case of a string with parabolic density 
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Fig. 4. Application of Theorem [2] to the Borg string with a = 10. The sets represent 
the ratio E'^/n^vr^ for the first five modes (n = 1,...,5) corresponding to the 
different iterations. 



p{x) 



[1 + ax) 



(71) 



where |a| < 2 implies that p{x) > for a; G (—1/2,1/2). The reason for 
studying this string is practical, since it allows a simple implementation of 
Theorem [3l 

We choose tiq^x) = {x+ 1/2) (L' — x) and apply eq. (!23l) to calculate the higher 
iterations, leaving L' as a shooting parameter. For example, after one iteration 
we find: 



ao + aix + a2X^ + a^x^ + a^x'^ 



where 



(72) 
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Fig. 5. Real part of the solutions to eq. ([73 
The condition rii{L) = 0, provides the quartic equation 



(a + 2){2L' - l)(a(4L'(L' + 1)(4L' + 1) - 5) + 10(4L'(L' + 2) - 1)) ={Q^) 

which has a fixed solution, corresponding to L' = 1/2 and three remaining 
solutions which depend on a. Although these last solutions can be calculated 
explicitly, their expressions are lengthy and not particularly enlightening: we 
prefer to give approximate expressions around a = 0: 



L' 



(74) 



L',^-{V5-2) (2a + 5) + 

^; 4a 5 3 ^ [ r- 

L'^ + - + a' 

^ 5 2a 4 L 

1 



a 



10 



2 + VS) (2a + 5) + O 



a 



In Fig. Owe draw the four solutions (just the real part). In order to be accept- 
able, the solutions must be real and —L < L' < L, i.e. they must fall in the 
region delimited by the two horizontal lines. The constant solution, L' = 1/2 is 
physically acceptable and corresponds to the fundamental mode of the string 
(i.e. to the Theorem [1]); for a > the acceptable solution is L[, while for a < 
the acceptable solution is L'^. 

Fixing a = 1 and working through order 20 we have been able to obtain 
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Fig. 6. |?72o(-^)| a function of L' . The arrows are placed in correspondence of the 
zeroes of 'r]2o{L) and locate approximate solutions of the Helmholtz equation. 

an explicit expression for ri2o{x): in Fig. [6] where we display |?72o(-^)| as a 
function of L', we see that there are 7 zeroes (see the arrows in the plot) 
which correspond to an equal number of approximate solutions of eq. ([1]). 
Notice that the rjjix) are not normalized and that we have used a logarithmic 
scale in the vertical axis (so that the zeroes of \rin{L) \ correspond to the spikes 
in the plot). 

In Table [3] we report the first seven eigenvalues of the Helmholtz equation 
([1]) corresponding to a parabolic density with a = 1 obtained using Theorem 
[3] with 20 iterations (third row). The second row reports the values of L' 
for which 7720 (-^) = 0; the fourth and fifth rows report the results obtained 
respectively with collocation with a grid corresponding to = 2000 and 
with the WKB formula. Notice that the first 5 eigenvalues obtained with the 
Theorem are in almost perfect agreement with the very precise results obtained 
with collocation; the last two eigenvalues on the other hand are rather poor. 

At this point we can make few observations: 

• The number of zeroes of f]j{L), j being the index of the iterations, grows 
with j: this means that by iterating a large number of times, one may 
obtain at once a large number of approximate solutions of equation ([T]); 
these solutions all correspond to the same rjj with a different value of L'; 

• No orthogonalization is needed, although it could still be used to improve 
the quality of the results; 

• Using a good ansatz may help reduce the number of iterations: 

• This approach requires that we can perform the integrals in eq. (1231) . which 
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7 


-0.246658591 


391.6830664 


482.3543658 


483.6106157 



Table 3 

First seven eigenvalues of the Helmholtz equation ([T]) corresponding to a parabolic 
density with a = 1 obtained using Theorem [3] with 20 iterations (third column). 
The second column reports the values of L' for which rj2o{L) = 0; the fourth and 
fifth columns report the results obtained respectively with collocation with a grid 
corresponding to iV = 2000 and with the WKB formula. 

become more and more involved as the iterations are increased; 



6 Conclusions 



In this paper we have discussed the problem of a string of variable density. 
Our results show that: 

• It is possible to build iteratively precise approximations to the fundamental 
and excited modes of the string, which converge to the exact ones as the 
number of iterations tends to infinity; 

• The results of perturbation theory to second order agree with those of the 
WKB method for the highly excited modes of the strings, and imply a 
functional dependence on the density which is different from the one in two 
or higher dimensions: in ref. [IT] we have obtained that Weyl's law for a 
d-dimensional cube of side 2L filled with density T,{xi, . . . , xj) is 

• A collocation method allows to obtain very precise numerical results for 
large number of modes; 



The author ackowledges support of Conacyt through the SNI fellowship. 



26 



References 



[I] M.G. Krein, Determination of the density of a nonhomogeneous symmetric cord 
by its frequency spectrum. Doklady Akad. Nauk SSSR (N.S.) 76, 345348. (1951) 

[2] M.G. Krein, On a generalization of investigations of Stieltjes. Doklady Akad. 
Nauk SSSR (N.S.) 87, 881884. (1952) 

[3] M.G. Krein, On inverse problems for a nonhomogeneous cord. Doklady Akad. 
Nauk SSSR (N.S.) 82, 669672. (1952b) 

[4] M.G. Krein, On some cases of effective determination of the density of an 
inhomogeneous cord from its spectral function. Doklady Akad. Nauk SSSR (N.S.) 
93, 617620. (1953) 

[5] M.G. Krein, On a method of effective solution of an inverse boundary problem. 
Doklady Akad. Nauk SSSR (N.S.) 94, 987990. (1954) 

[6] R. Beals, D.H. Sattinger and J. Szmigielski, The string density problem and the 
CamassaHolm equation, Phil. Trans. R. Soc. A 365, 2299-2312 (2007) 

[7] J.W.S. Rayleigh, The Theory of Sound vol 1 (New York: Dover) (1945) 

[8] G. Borg, Eine Umkehrung der Sturm-Liouvilleschen Eigenwertaufgabe. Acta 
Math. 78, 1-96. (doi:10.1007/BF02421600) (1946) 

[9] C.O. Horgan and A.M. Chan, Vibration of inhomogeneous strings, rods and 
membranes. Journal of Sound and Vibration 225, 503-513 (1999) 

[10] H.RW. Gottlieb, Isospectral strings. Inverse Problems 18 971-978 (2002) 

[II] L.P. Fulcher, Study of the eigenvalues of a nonuniform string, Am.J.Phys.53 
730-5 (1985) 

[12] C.M. Bender and S.A. Orszag, Advanced mathematical methods for scientists 
and engineers, McGraw-Hill (1978) 

[13] F.S. Crawford, Eigenvalues of a nonuniform string, Am.J.Phys. 55, 168-169 
(1987) 

[14] C. Castro and E. Zuazua, High frequency asymptotic analysis of a string with 
rapidly oscillating density, European Journal of Applied Mathematics 11, 595-622 
(2000) 

[15] C. Castro and E. Zuazua, Low Frequency Asymptotic Analysis of a String with 
Rapidly Oscillating Density, SIAM J. Appl. Math. 60, 1205-1233 (2000) 

[16] G. Rawitscher and J.Liss, The vibrating inhomogeneous string: a topic for a 
course in Computational Physics, l arXi v:1006.1913vl (2010) 

[17] P. Amore, Can one hear the density of a drum? Weyl's law for inhomogeneous 
media. larXiv:0912.1402v l [math-ph] (2009) 



27 



[18] P. Amore, Solving the Helmholtz equation for membranes of arbitrary shape: 
numerical results, Journal of Physics A 41, 265206 (2008) 

[19] P. Amore, Spectroscopy of drums and quantum billiards: perturbative and non 
perturbative results, J. Math. Phys. 51 (2010) 

[20] C. Alvarado and P. Amore, Spectroscopy of annular drums and quantum rings: 
perturbative and non perturbative results, arXiv:1004.2407 vl [quant-ph] (2010) 

[21] P. Amore et al., Variational collocation on finite intervals. Journal of Physics A 
40, 13047 (2007) 

[22] P. Amore et al.. Collocation on uniform grids. Journal of Physics A 42, 115302 
(2009) 



28 



